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Abstract 

We study synchronization phenomena in scale-free networks of asymmetrically coupled dynamical 
systems featuring degree-degree correlation in the topology of the connection wiring. We show that 
when the interaction is dominant from the high-degree to the low-degree (from the low-degree to the high 
degree) nodes, synchronization is enhanced for assortative (disassortative) degree-degree correlations. 
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Recent years have witnessed increasing interest in the study of the structure of complex networks. Networks 
are present ubiquitously in the real world, including biological systems, genetic chains, social relationships and 
artificial and engineering architectures A main quantity that characterizes the structure and dynamics 
of those networks is the so called degree distribution -P(fc), the probability that a randomly chosen node 
within the network has degree k, i.e. there are k connections incident in it. In most cases, this distribution 
has been found to obey a power law P{k) k^'^, with 7 usually ranging between 2 and 3. The slow 
decay of the power-law distribution, when the exponent is within the observed values, is responsible for the 
presence of few high-degree nodes, usually termed as hubs, which play a leading role in characterizing the 
network dynamical behavior. Thus the network structure (termed in this case as scale free) displays high 
heterogeneity, with the presence of a few high-degree vertices and many low-degree ones. 

Beyond the specific interest in the study of the main topological properties of such networks, there is an 
ongoing research effort to understand how these statistical features affect the dynamical processes taking 
place over them, as e.g. epidemic spreading [531131 [Hj, communication models [531 [TUl [HI [I] and spin systems 
[31 [HI HI]- Among these, a relevant phenomenon is the emergence of collective synchronized motion. This 
corresponds to the case of several dynamical systems, starting from different initial conditions appropriately 
coupled with a network topology, stabilizing in a common evolution. 

In this framework, we consider the case of networks of coupled identical systems, that can be conveniently 
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described by the following set of ordinary differential equations: 

^ = fix,) -af^C^M^j), z - 1, 2, ...N, (1) 

where Xi G M™ is the m-dimensional vector describing the state of the i'^ node, f{xi) : R™ — > R™ governs 
the dynamics of the z*'' node, h{x) : M™ R™ is a vectorial output function, ct is a parameter ruling the 
coupling strength, and N is the number of nodes of the network. The topological information on the graph 
is contained in the connection matrix C. Specifically C = {Cij} is defined as follows: 

Cij = -AijWij, i ^ j 
j 

where A is the irreducible adjacency matrix associated to the network topology, that is Aij = Aji = 1, if 
there is a connection between i and j, otherwise. W is a positive definite matrix containing the information 
about the weights over the network connections, with Wy giving a measure of the strength of the interaction 
from node i to node j. 

Note that ^ and ([2]) encompass a general and large class of networks, allowing for non-linearity in the 
functions / and h and leaving complete freedom for the choice of the weights over the network connections. 
The only constraint is represented by the condition Cu = — X^j^^ i ^ij ^ which is necessary in order to guarantee 
that the connection matrix C is zero row-sum. 

When the connection matrix C is symmetric, the network propensity to synchronization (the so-called net- 
work synchronizability) can be studied rigorously by means of the Master Stability Function (MSF) approach 
|25j . Namely, the MSF assesses the necessary conditions for the transverse stability of the synchronization 
manifold xi = X2 — ■■■ = xn ^ Xg, whose existence as an invariant set for the system in ([1]) is warranted by 
the zero row-sum property of C 

Indeed, let Sxi{t) — Xi{t)—Xs{t) be the deviation of the i^^ vector state from the synchronization manifold, 
and consider the m x N column vectors X — [xi, X2, ■ ■ ■ , xn)'^ and (5X = (Sxi, . . . , Sxm)'^ ■ 

Then, one has 



6± = [lN®Jf{xs)-cTC®Jh{xs)]51^, (3) 

where J f and Jh represent the Jacobian operators of the system dynamics and coupling. Following [25], 
Eq. ([3]) can be block-diagonalized in the form: 

'^ = [J.f{xs)-<y\jJh{xs)]m. j = 1,2, ...,7V, (4) 

where {Aj} = {A^ + zA*} is the set of eigenvalues of the coupling matrix C (supposed diagonalizable) , 
ordered in such a way that A^^ < X^--- < A^. Notice that Ai is identically equal to 0, due to the zero row- 
sum condition, and the corresponding eigenvector is the unitary vector pointing in a direction parallel to the 
synchronization manifold. Therefore, the modes of the system along this direction are not relevant, since one 
has only to care for whether the system will synchronize or not, independent of where this happens in phase 
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space. Thus, the stabihty of the synchronization manifold should be checked with respect to perturbations 
lying in phase space along directions orthogonal to it (those associated to the eigenvalues A2, Ajv). 

In order to study the stability of the synchronization manifold as a function of the eigenvalues of the 
matrix one introduces the following parametric equation in a: 

^[J.f{xs) -aJh{xs)]r], (5) 

where a takes the place of aXj in eq. Note that Jf and Jh do not depend on the network topology. 
Moreover, for simplicity, one can rewrite ^ as: 

§ = K{aH (6) 

where K{a) is the evolution kernel, as a function of a. The MSF associates to each value of the parameter 
a in the complex plane, the maximum (conditional) Lyapunov exponent of the system calculated along 
trajectories of ([6]). 

Now, in most cases (i.e. for most of the forms of the functions / and h), the MSF is negative in a bounded 
region of the complex plane (termed as the stability region S) . The immediate consequence of this is that 
the network synchronizability is related to the distribution of the eigenvalues of the coupling matrix C 
Specifically, the stability of the synchronization manifold requires as a necessary condition that all the aXj , 
j = 2, ...,N belong to S. 

As proposed in [12,, a sufficient condition to verify the above requirement is checking that the rectangle 
of the complex plane, TZ: [crA2,crA^] x [crminj A*, cr maxj A*] is fully contained in S (observe also that the 
spectrum of C is symmetric with respect to the real axis and thus M — max^ A* = — minj A* ) . Moreover 
this condition can be simplified with respect to different forms of the functions / and h into the general 
requirement that both the eigenratio defined as i? = ^^n/-^2 Af, the maximum imaginary part of the 
spectrum of C, are as small as possible. In so doing, different network topologies become directly comparable 
in terms of their synchronizability, by simply looking at the indices R and M . 

Moreover, we wish to emphasize that there is no need for introducing an ordering between the two 
parameters R and M, defined above; specifically, when they are both reduced, this leads to the general 
result that the range of values of cr, say I^^, for which TZ is contained in S (and therefore the network 
synchronizability) increases. Note that this result is independent from the specific forms of / and h, i.e. of 
the particular MSF considered. 

It is worth noting that the block-diagonalization in (|4]) can be properly applied only to the case of sym- 
metric and asymmetric diagonalizable matrices (in the particular case of a real spectrum, as e.g. when C is 
symmetric, we obtain the following simple relationship: the lower is R, the larger is the interval I^). More- 
over, as explained in [211 [26], the same definition of network synchronizability (in terms of both R and Af) 
can be extended also to the case of non-diagonalizable matrices (see e.g. [HI [TBI [U [12] ) , when the condition 
is satisfied that the network embeds at least an oriented spanning tree. The case of generic directed networks 
has been studied in [5]. 

Thus in what follows we shall seek to characterize how the structural properties of complex networks, in 
terms of both the forms of the adjacency matrix A and the weights matrix W, may affect (and eventually 
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improve or hinder) their synchronizabihty. 

The first contribution in this sense was provided in [22] , where unweighed scale free networks were compared 
in terms of their synchronizabihty as varying the degree distribution exponent 7. 

Note that in the case of unweighed topologies, W is assumed to be the unitary matrix and thus £ can be 
simply expressed as £.y = — Aj for i ^ j and La = J2j -^iji * — l,2,...,Af. 

It can be shown that as the exponent 7 is increased, the network topology becomes more homogeneous and 
the average distance within pairs of connected nodes along the geodesic increases. For instance, real-world 
networks are characterized by low values of 7 in that they are affected by high heterogeneity in the degree 
distribution and low distances between pairs of connected sites (the so-called small-world effect). 

Thus one would expect such networks to be characterized by high synchronizabihty, in order also to justify 
their onset in nature as being the result of a self-organized process. Instead, the analysis in [52] showed that 
these are typically characterized by lower synchronizabihty than their homogeneous counterparts, giving rise 
to an unexpected phenomenon, which is termed in the literature |16j as the paradox of heterogeneity. 

This first result was overturned in [15j where an appropriate choice of the weights was considered over the 
network links of the form: 

W,,=k-^, i^j, (7) 

where ki is the degree of node i, i.e. ki — Aij and /? is a variable parameter. 

Now, it is worth noting that in the real world, weights associated to the network links depend on an huge 
amount of variables; furthermore they are typically also variable in time. Thus, trying to estimate the effects 
of such a complex behavior, may be a difficult task. On the other hand, the approach proposed in |15| is 
very convenient in that it is aimed at evaluating the effects of variable weights, depending essentially on 
some network topological features. Namely, by varying /3, it is possible to tune the strength of the coupling 
among the network nodes according to their degree: f3 — Q corresponds to the case of unweighed topologies, 
positive (negative) values of (3 indicate that the strengths of the couplings acting on each vertex decreases 
(increases) with its degree. In |15| it was shown that the network synchronizabihty is optimized at /3 = 1, 
whatever the form of the degree distribution. Specifically, it was claimed that in this optimal regime, the 
eigenratio is independent from the networks topological properties, other than the average degree 

(note that for a coupling of the form in ([7]), the parameter M = [T5]). 

The case of /? = 1 is of particular relevance also in another aspect. In fact, in such a case, Cu — 1 Vi and 
J2j^i^ij ^ ^1 "^^^ Gerschgorin circle theorem [T^] states that the set of the eigenvalues {Xj = + jXj}, 
belongs to the union of the circles (Ci) centered on the main diagonal values of C in the complex plane 
and with radius equal to the sum of the absolute values of the other elements in the corresponding rows, 
{{Xi} C C\iCi[Cii,J2j^i Thus in the case of /? = 1, all the d coincide with the circle of radius 1, 

centered at 1 on the real axis, inside which are constrained to lie all the Xj, j = 1, 2, N ||]. Note that 
from this condition, we can draw the two following ones: (i) = < A2 < ... < A^ < 2, (ii) A* £ [—1, 1], 
j = 1,2, ...,7V. 

It is worth noticing here that for any given choice of the weights over the network connections , the 
same requirements as above are satisfied by replacing ^ with 
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A, = 1, i=J- (8) 

We wish to emphasize that this can be convenient in comparing networks characterized by different 
topological features in terms of their synchronizability. In fact, although the condition on the minimization 
of corresponds to an increase of the interval of the values of a for which the synchronization manifold 

is stable, it does not give any information about the values of a themselves. One possible consequence of 
this could be e.g., that although a network is more synchronizable than an other (in terms of their respective 
eigenratios), it experiences synchronization only for higher values of the coupling a. Nonetheless this problem 
can be overcome by considering an appropriate coupling as in ([5]) , which guarantees that A2 is lower-bounded 
by and Aat is upper-bounded by 2 and thus ensures the spectra we are evaluating are characterized by the 
same scale. 

By following this approach, in ^ it was shown that an even better result in terms of network synchroniz- 
ability may be achieved by taking weights proportional to a specific topological quantity known as load [lOj 
or betweenness centrality [17j . At the same time in |,12 , a particularly suitable way was proposed of choosing 
the weights over the network connections, according to the presence of a given ordering among the network 
nodes. Specifically the network nodes were ordered according to a given scalar property, say p, in such a 
way, i.e. that i < j if p{i) < p{j)- Then the following general way of choosing the weights was proposed: 



W., = ^ 7^ , (9) 

where Qij = {Qij = ^^^) if * < J (* > j)i with —1 < g < 1 being a parameter governing the coupling 
asymmetry within the network. Note that, by varying it is possible to study the effects of variable 
asymmetric weights over the network edges. Specifically, in what follows we will consider a very natural 
way of choosing the weights based on a degree ordering (also termed as age ordering in [12]), i.e. such that 
p{i) = ki. Namely, q = Q corresponds to the case of symmetric coupling, which is equivalent to the choice 
of /3 = 1 in eq. ([7]). The limit oi q ^ 1 [q ^ —1) corresponds to the case of a direct network in which 
only the connections from high-degree to low-degree vertices (from low-degree to high-degree) are present. 
Under this condition, it was shown in [12] . that in the case of scale-free networks, synchronizability can be 
improved for decreasing values of g, i.e. when the strength of the coupling is dominating in the direction 
from the high-degree to the low-degree vertices. 

Sofar, we have considered the effects of varying the weights over the network connections on the network 
synchronizability (in terms of the matrix W), while the network topology (in terms of the matrix A) has 
been considered as fixed. Nonetheless, it was shown in [22l[6], that even in the case of unweighed topologies, 
synchronizability may be strongly affected by varying some particular network topological properties. 

The most important one is, as already stated before, the degree distribution, which affects in a non- 
negligible way almost all the dynamical processes taking place over networks 2 . On the other hand, many 
important properties have been discovered to characterize in more detail the network structure, other than 
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the degree distribution. These are mostly due to particular forms of correlation or mixing among the network 
vertices [20] . 

One form of mixing is the correlation among pairs of linked nodes according to some properties at the 
nodes. A very simple case is degree correlation [18) . in which vertices choose their neighbors according to 
their respective degrees. Non-trivial forms of degree correlation have been experimentally detected in many 
real-world networks, with social networks being typically characterized by assortative mixing (which is the 
case when vertices are more likely to connect to other vertices with approximately the same degree) and 
technological and biological networks, by disassortative mixing (that takes place when connections are more 
frequent between vertices of different degrees). In [18] this property has been conveniently measured by 
means of a single normalized index, the Pearson statistic r defined as follows: 



where qk is the probability that a randomly chosen edge is connected to a node having degree fc; aq is 
the standard deviation of the distribution qk and Ckk' represents the probability that two vertices at the 
endpoints of a generic edge have degrees k and k' respectively. Positive values of r indicate assortative 
mixing, while negative values characterize disassortative networks. 

In [B] a strategy was used, similar to the one presented in [18| , which allows to vary the degree correlation (in 
terms of variable values of r) , while keeping fixed the degree distribution and it was shown that disassortative 
mixing, i.e. the tendency of high-degree nodes to establish connections with low-degree ones (and viceversa), 
is indeed a desirable network property in terms of its synchronizability. The same effect was checked in [7] 
to be persistent under a normalized form of the weighing as that proposed in |15| . with (3=1. 

Here, we shall seek to characterize how the combination of the effects of variable degree correlation in 
terms of the structural configuration of the network, and variable asymmetry over the network weights may 
indeed affect the network synchronizability. 

In order to construct networks, characterized by a given degree distribution (in terms of the exponent 7), 
we consider the network building model proposed in 13J, which is commonly termed as the configuration 
model. Specifically this consists in the following steps: 

1. Start from a given form of the degree distribution P{k). 

2. Assign to each of N vertices a target degree k > kmin drawn from the distribution P{k). 

3. Randomly create connections between pairs of vertices, avoiding multiple and self-connections, until 
all vertices have reached their target degree. 

Then in order to reproduce a desired level of degree correlation, following |20j . we propose exchanges of 
links between pairs of connected nodes, until the desired value of the observable r is reached. 

In what follows we have evaluated the effects on the network synchronizability of both variable values 
of r and q. The main results are shown in Figs. [T] and [21 for scale free networks with degree distribution 
exponent 7 = {2,3}. 




1 k,k' 



(10) 
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Figure 2: Maximum imaginary part of the spectrum of £, M, as function of q for scale free networks 
characterized by variable degree correlation: r = [—0.3,-0.15,0,0.15,0.3]. Left panel: 'y — 2, N ~ 1000. 
Right panel: 7 = 3, iV = 1000. 

The effects on the eigenratio R have been reported in Fig. [T] Note that, as already shown in [12] for 
uncorrelated networks, R is an increasing function of q for all the values of r considered. This depends on 
the advantage, in terms of the network synchronizability, of having asymmetric interactions directed from 
high-degree nodes to low-degree ones. Moreover this is particularly evident for low values of q in the case 
of assortative networks. Thus the combination realized when q ^ —1 and the network topology is strongly 
assortative (high r), represents the optimum configuration for the minimization of R. 

This is also consistent with the qualitative claim presented in [12 that the optimal network configuration 
is obtained, when both a dominant interaction from high-degree nodes to low-degree ones and a structure 
of interconnected hubs are present. 

On the other hand, when 5 « 0, we observe a completely different picture with disassortative networks 
being characterized by better synchronizability properties (as already shown in :7J). The same behavior is 
also confirmed in the case of positive values of g, where the already better performance of disassortative 
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Figure 3: Spectrum of networks characterized by different degree correlation properties and degree of asym- 
metry. Each picture represents A* vs A^, j = 1, 2, N. Green points are used for the spectrum of networks 
characterized by negative degree correlation (r = —0.3), black points are used for the spectrum of networks 
characterized by absence of correlation (r = 0), red points for the spectrum of networks characterized by 
positive degree correlation (r = +0.3): {&) j — 2, q ~ 0.8; (b) j — 2, q = —0.8; {c) j = 3, q — 0.8; (d) 7 = 3, 
q = -0.8. 

networks is further enhanced through asymmetric coupling. 

Thus the onset of two separate regimes emerges as varying q: the first one, in the case where there is a 
dominant interaction of the high-degree nodes on the low-degree ones, with assortative mixing representing 
a desirable property in terms of the network synchronizability; the second one, with symmetric coupling or 
directed from low-degree nodes to high-degree ones, where disassortative mixing represents the best choice 
in order to enhance the network synchronizability. 

The behavior of the maximum imaginary part of the spectrum, M, as varying q, is show in Fig. [2l As 
expected, this is characterized by a minimum at q = (for which the coupling is symmetric) and increases 
as increasing the asymmetry over the network links. Note that for all values of g / 0, assortative (positively 
correlated) networks are characterized by higher values of M when compared to their uncorrelated and 
disassortative counterparts. 

A more exhaustive picture is provided in Fig. O Specifically here two configurations characterized by 
strong asymmetry {q = 0.8 and q — —0.8) are compared for scale free networks with exponents 7 = {2,3}, 
in terms of their spectra: AJ vs A*, for j = 1, 2, .., A^. Networks characterized by different degree correlation 
properties are depicted with different colors: red for assortative mixing, green for disassortative, and black for 
no-correlation. Remind that the condition for increasing the network synchronizability is that the spectrum 
{Xj} for j = 2, .., N is contained in a region of the complex plane as small as possible. 

This condition is clearly achieved for negative values of q (compare respectively Figs. [Sja) with [31(b) 
and Figs. [31(c) with[3l^d)). Note that in the case oi q — 0.8, the advantage of introducing negative degree 
correlation is undoubted, as both R and M are reduced and the green area is clearly smaller than the red 
one in Figs. [31(a) and [31(c). On the contrary, in the optimal regime where q = —0.8, the reduction of the 
range of the real values of the spectrum A^' , j = 2, 3, A^ is counterbalanced by an increase in the range of 
the imaginary parts A*, j = 2, 3, A^, so that a more intricate phenomenology emerges. 
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Figure 4: Synchronization error {E) as a function of the couphng strength a for a scale- free network of 1,000 
Rossler oscillators (configuration model, 7 = 2). Evolution equations and parameters are specified in the text. 
Each point corresponds to an ensemble average over a set of 10 independent realizations. The asymmetry 
and correlation parameters are q — 0.8, r = 0.3 (diamonds), q — 0.8, r — —0.3 (circles), q = — 0.8, r — 0.3 
(triangles) and q = — 0.8, r = —0.3 (squares). The numerical simulations support the qualitative scenario of 
the MSF: for positive (negative) values of q, synchronization is enhanced in disassortatively (assortatively) 
correlated graphs. 

Specifically, in such a case it becomes necessary to look at the specific form of the MSF under evaluation, 
in order to state the benefit of introducing degree correlation among the network connections. This should 
be taken into account, also with respect to variable forms of the MSF, when different choices of the functions 
/ and h are taken into consideration. 

In order to confirm this qualitative description with an example, in the following we numerically simulate 
a network of 1,000 coupled chaotic Rossler oscillators. The dynamics of the system is ruled by Eq. [1] with 
m — 3, X = (xi, X2, X3), f{x) = [—X2 — X3, xi + 0.165x2, 0.2 + X3{xi — 10)] and h{x) — [xi, 0, 0]. 

The appearance of the synchronous state can here be monitored by looking at the vanishing of the time 
average (over a window T) of the synchronization error (E) = "^jyi It^^ I l^i I \dt' . In the present 

case, we adopt as vector norm = \xi \ + \x2 \ + \x3\. 

Fig. m reports (E) vs. a for different networks topologies and degree correlations. It is easy to observe 
that the numerical simulations fully support the qualitative scenario of the MSF. Indeed, by comparing the 
synchronization effects in the cases q = 0.8, r = 0.3 (diamonds), q — 0.8, r = —0.3 (circles), q = —0.8, r = 0.3 
(triangles) and q = — 0.8,r = —0.3 (squares), one immediately realizes that for positive (negative) values of 
q, synchronization is enhanced in disassortatively (assortatively) correlated graphs. 

In conclusion, in this paper we have regarded the network topology as the set of active connections 
present in the network, that is pairs {ij}, i,j = 1...N, such that Aij = 1. We analyzed the effects of 
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varying the weights over the network connections (by making the parameter q change) in order to investigate 
different levels of asymmetry. Moreover, we studied how variable degree correlation, while keeping the 
degree distribution unchanged, can favor or hinder the overall synchronization. While negative degree 
correlation was found to be a convenient network structural property in terms of synchronizability in the 
case of symmetric coupling {q = 0) and of positive values of g , a more intricate phenomenology was observed 
in the optimal situation, i.e. for low values of q. Specifically, in such a case an interesting balance between 
the width of the range of the real and imaginary parts of the eigenvalues has been observed. 

We wish to emphasize that this issue could be of relevance in those applications where, because of different 
circumstances, synchronization may be or not be a desirable effect. 
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